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The non-linear hydrodynamic equations for axion/scalar field dark matter (DM) in the non- 
relativistic Madelung-Shcrodinger form are derived in a simple manner, including the effects of 
universal expansion and Hubble drag. The hydrodynamic equations are used to investigate the 
relative velocity between axion DM and baryons, and the moving-background perturbation theory 
(MBPT) derived. Axions massive enough to be all of the DM do not affect the coherence length of 
the relative velocity, but the MBPT equations are modified by the inclusion of the axion effective 
sound speed. These MBPT equations are necessary for accurately modelling the effects of axion DM 
on the formation of the first cosmic structures, and suggest that the 21cm power spectrum could 
improve constraints on axion mass by up to four orders of magnitude with respect to the current 
best constraints. A further application of these results uses the “quantum force” analogy to model 
scalar field gradient energy in a smoothed-particle hydrodynamics model of axion DM. Such a model 
can treat axion DM in the non-linear regime and could be incorporated into existing N-body codes. 


I. INTRODUCTION 

The particle nature of the dark matter (DM) is 
unknown, yet cosmological and astrophysical probes 
provide a wealth of information about the length 
and time scales over which it must form structure, 
providing bounds on the DM particle mass. In the 
case of thermal DM, structure formation excludes 
hot DM and limits warm (W)DM to have mass 
mw ^ C>(few keV), e.g. Refs. pQ [5]. In the case 
of non-thermal DM, for example ultra-light axions 
(ULAs) or other scalar fields, e.g. Refs. [MS], the 
limit on the mass is rria ^ O(few) x 10“^^ eV, e.g. 
Refs. [HHH]- Forthcoming astrophysical data, for ex¬ 
ample from the 21cm power spectrum, could extend 
these bounds substantially and possibly find evi¬ 
dence pointing to the particle nature of DM. Util¬ 
ising this data requires models for DM on small 
scales, handling non-linear effects. WDM and self¬ 
interacting DM are fairly well studied in this regard, 
e.g. Refs. [5HTT]. but much less is known about how 
to model ULAs and other scalar helds. The most ad¬ 
vanced simulations in this held are those of Schive 
et al [12] for studying dwarf galaxy sized objects 
and density cores. Alternative simulation models are 
necessary, both to conhrm the accuracy of Ref. HH 
simulations and to study new effects on different 
scales. 
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In this note I present a formalism, the non¬ 
linear hydrodynamic equations (NLHEs) in the non- 
relativistic limit. This is a useful model to com¬ 
pute the effects of ULAs on the 21cm power spec¬ 
trum. The cut-off of power at small scales, k > 
0{1) Mpc“^, induced by ULAs satisfying rria > 
10“^^ eV, via the effect on star formation and the 
baryon-DM relative velocity nMm, has a knock-on 
effect at large scales that could distinguish ULAs 
from CDM. The NLHEs also provide a possible 
method to model ULAs using smoothed-particle 
hydrodynamics (SPH), which is distinct from the 
method of Ref. m who use the Widrow-Kaiser |19] 
approach. The SPH method can be more easily in¬ 
corporated into existing N-body codes, such as GAD¬ 
GET [^, which already contain SPH modules to deal 
with baryon pressure. The method has been dis¬ 
cussed recently in Ref. |21j , and draws on ideas from 
Bohmian mechanics and quantum trajectories (e.g. 
Refs, [53]) to compute the anomalous pressure 
due to the scalar held gradient energy. I present 
the equations, and discuss various issues and links 
to other studies; possible algorithms are discussed 
in Refs. mills], and will be the subject of future 
discussion. 


All linear power spectra are computed using AX- 
ionCAMB [8], a modihed version of the publicly 
available CAMB code [H]. I work in units where 
h=c=l. 



II. HYDRODYNAMIC EQUATIONS 


The action governing an axion (or general scalar 
field([^ of mass rria in general relativity is 


S' = y 




( 1 ) 


Assuming the vanishing of anisotropic stress, the line 
element in Newtonian gauge in an expanding uni¬ 
verse is 


ds^ = -{l + 2V)dt^+a^{l-2V)dx^ , (2) 


the equation for the WKB amplitude '0 (and also 
independently for 0*): 

itp — 3iiS0/2 -I- (2TOaa^)“^V^0 — maVtp = 0 . (6) 

Eq. ([^ is the Schrodinger equation on an ex¬ 
panding spacetime. The Madelung form (see below) 
of this equation has been presented before in e.g. 
Refs. 123 HI] by writing the standard Schroodinger 
equation in the comoving frame. Expressing the 
equations in this frame is not the norm in studies of 
the Schrodinger-Possion system [UlIlHHj, where 
expansion effects are accounted for by rescaling the 
wavefunction by the background solution. 

The axion energy density, is found from the 
energy momentum tensor 


where a is the cosmic scale factor and V is the New¬ 
tonian potential. The Hubble rate H = d/a. The 
scalar held 0 obeys the Klein-Gordon (KG) equation 

(n-r7i2)0 = O. (3) 

To hrst order in U ~ (weak-held limit, where 
enr is the perturbative parameter for relativistic ef¬ 
fects) the D’Alembertian is 


^(5“'^5a00/30 + m^0^), (7) 

with r°Q = -pa giving 

Pa = ^[(1 - 2E)0^ -k -k a"^(l -k 2V)d'’(j)di(j )\. 

( 8 ) 

to hrst order in V. Taking the same limits as before 
we hnd the leading order piece is 

p, = |0p+O(e). (9) 


□ = -{l-2V){d/ + iHdt)+a-‘^{l + 2V)V‘^-Wdt . 

(4) 

In the limit that H/ma ~ ewkb we can use the 
WKB approximation to solve for 0, giving 

0 = (m,v^)-i(0e-™“‘+0*0™“‘). (5) 

Taking WKB to hrst order implies 0/ma0 ewkb- 
I also take the non-relativistic limit by letting 
the dispersion relation for wavenumber fc be w = 
\/ml + -k kml + C>(e^j^), i.e. k/rua ~ enr 

|26| . Performing a double expansion to 0(e^j^ wkb)’ 
the non-relativistic limit of the KG equation gives 


Cosmological perturbation theory for scalar helds 
normally makes a background-huctuation split at 
the level of the held, 0 = 0-k<50 (e.g. Refs. [23|30]). 
This, however, does not preserve the canonical form 
of the huid equations at non-linear order in den¬ 
sity huctuations. Therefore, I make the background- 
huctuation split at the level of the density, as is usual 
for CDM and baryons (e.g. Ref. [?T]1. 

To do this, express the wavefunction, 0, in polar 
co-ordinates as 0 = Re*'® (known as Madelung form) 
and transform Eq. ([^ into equations of motion for 

Pa = R^, Va = imaa)~^VS, ( 10 ) 

where this defines the axion fluid velocity, Va- Sep¬ 
arating real and imaginary parts of Eq. ^ gives 


^ I do not include self-interactions, as they are model- 
dependent. One can show that they are sub-dominant to 
gravity on linear and non-relativistic scales for an axion 
with the canonical cosine potential. Scalar field conden¬ 
sates with attractive and repulsive self interactions are dis¬ 


cussed in e.g. Ref. m- Axions have an attractive quartic 
self-interaction, and thus no Thomas-Fermi limit. Other 
aspects of the self-interaction will be discussed in future 
work. 
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( 11 ) 


Pa + SHpa+ — {paVa) = 0 , 


Hv ■ 



1 - ^ i 

U'a • — 

]Va = - 

V a ) 

' a \ 


1 

2maO? ,/pa 


( 12 ) 


The background-fluctuation split can now be carried 
out as usual by writing Pa = Pa{^ + 5a)- This gives 

Pa T ‘5Hpa — 0 ^ Pa — HaPcrit (a/uo) 5 (1^) 


i.e. the axion background energy density behaves, 
in the WKB approximation, as a w = 0 barotropic 
fluid (true for an oscillating field with a harmonic 
potential minimum [32)1. The fluctuations obey 


5a + a-Va • V5a = -a-\l + 5a)V • Va , (14) 

^a + a-^ (va • V) Ua = -a-^V{V + Q)-Hv, (15) 
1 


Q = - 


2mla^ Vl +5a 


(16) 


where I have defined the “quantum potential” Q that 
accounts for the scalar-field gradient energy. When 
Q —>■ 0, scalar fields have a Jeans instability [53]. 

A background-fluctuation split on the Einstein 
equations gives 


Pi, VV = 47rGa^^p,Ji, (17) 


where the second equation applies in the sub¬ 
horizon, H/k ^ cnr, limit (see, e.g.. Ref. [34]), and 
the sums extend over all species, i. 

Eqs. (13-{^) are the complete NLHEs in the non- 
relativistic limit for axion/scalar-field DM. These 
equations differ from those for CDM only by the in¬ 
clusion of the Q term in the Euler equation, Eq. ®. 
They are valid in the non-relativistic limit even for 
non-linear density and axion field perturbationsj^ 
They thus illustrate the equivalence between axion 
DM and CDM on large scales (where k/rua 1), to 


^ The limits taken are valid after axions begin oscillating 
(H < rria), on sub-horizon but super-Compton length 
scales (H < k < rria), and for weak-fields (|i/| < 1). The 
first two of these are generally valid during structure for¬ 
mation in the matter dominated era down to sub-kiloparsec 
scales for the allowed axion DM particle masses. The weak- 
field limit can be violated, of course, via the instability to 
black hole formation. 


all orders in Newtonian perturbation theory. They 
also provide the correct setting to apply all the tools 
of standard perturbation theory m and modern de¬ 
velopments such as the effective field theory of large 
scale structures [5S| to axion DM. 

III. APPLICATIONS 

A. Relative velocity of axion DM and baryons 

The relative velocity between CDM and baryons 
can have an observable effect by suppressing star 
formation in the first structures at high-z UMI]. 
The inclusion of the quantum pressure, Q, in the 
Euler equation for axion DM demands a separate 
treatment. The relative velocity between axion DM 
and baryons at recombination, Zrec ~ 1020, can be 
computed in first order cosmological perturbation 
theory using a Boltzmann code. 

The variance of the relative velocity, Vba = vt — Va, 
is 



which defines the relative velocity power spectrum, 
^vba^ from the velocity divergence, 9i = a“^V • Vi- 
The primordial curvature power spectrum is = 

( \ rig —1 

^ j with amplitude, Ag, and tilt, rig, which 

are both well measured by Planck [36j : the pivot 
scale, fco, is conventional. Fig. shows for var¬ 
ious axion masses. The lightest axions suppress the 
relative velocity power at large wavenumbers. On 
the scales shown, TOq = 10“^^ eV is indistinguish¬ 
able from CDM. 

Relative velocities are gauge invariant, but the ve¬ 
locity of an individual species depends on the choice 
of gauge. For example, the CDM velocity divergence 
vanishes in synchronous gauge, 0c(syn.) = 0. Never¬ 
theless it is instructive to look at the single-species 
velocity power. Fig.[^shows the velocity power com¬ 
paring CDM and an axion with nria = 10“^^ eV in 
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FIG. 1. Relative velocity power spectrum for axions and 
baryons (Eq. |18[ ) at Zrec = 1020 for various axion masses. 
On the scales shown the rua = eV spectrum is 

indistinguishable from CDM. 



FIG. 2. Velocity power spectra for axions, rria = 
10“^^ eV, and CDM in synchronous and conformal New¬ 
tonian gauges. Note that the synchronous gauge has 
Uc(syn) = 0. The gauge terms decay on sub-horizon 
scales. 


synchronous and Newtonian gauges [31]. We ob¬ 
serve that, in the Newtonian gauge axion DM has 
suppressed velocity power with respect to CDM for 
wavenumbers k > kj^a, where kj^a = 1.6a\/Hma is 
the axion Jeans scale. In the synchronous gauge, 
where Va is the relative velocity between CDM and 
axions, axions have velocity power for k > kja- For 
TTia < 10“^"^ eV the coherence scale of the axion- 
CDM relative velocity could be large and may serve 
as an additional probe of the composition of the DM, 
complementary to the CMB constraints of Ref. |5] 
(see Refs. [371 [3H] for the case of neutrino-CDM rel¬ 
ative velocity). 

The baryon-CDM relative velocity, Vbc, is coher¬ 
ent on scales kcoh. ^ 1 Mpc“^. Large scale struc¬ 
ture and the CMB [T] |H] constrain ULAs to have 
TTia > 10“^^ eV if they are to be all of the DM. This 
implies that kj^a > ^coh. and axion DM has the same 
large scale velocity relative to the baryons. Thus, 
moving background perturbation theory (MBPT), 
which treats the large scale relative velocity non- 
perturbatively, can be carried out for axion DM just 
as for CDM [T3] for cell sizes of order a few Mpcs. 

MBPT works by writing the velocities as 

Vi{x,t) = + Ui[x,t), (19) 

then moving to the frame where 

= 0 and = —i^a®^ • (20) 

and performing first order perturbation theory on 
the variables 5i,Ui. Going to Fourier space and ex¬ 
panding to first order in 5a the quantum pressure, 
Q, is 

+ ( 21 ) 

The equations of motion for the coupled axion- 
baryon system in MBPT are therefore 


L=ia-\v^ba\t)-k)5a-ea, ( 22 ) 

Oa = ■ k)0a " 3H^{na{t)5a + VLb{t)5b)/2 - , (23) 

4 = -Ob , (24) 

Ob = -iH^{na(t)5a + nb{t)Sb)/2 - 2Heb + clbk^Stla^ , (25) 
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FIG. 3. Total matter power at 2 = 40 with (solid lines) 
and without (dashed lines) the Vbc effect, computed for 
various axion masses. For ma = 10“^* eV axion sound 
speed and vtc have effects on comparable scales. 


where Cg h is the baryon sound speed. In these equa¬ 
tions we identify the axion effective sound speed: 


c 


2 

s,a 




(26) 


which can be seen here as the manifestation of quan¬ 
tum pressure for linear density perturbations. 

The axion effective sound speed is the only term 
that distinguishes axion DM from CDM on small 
scales, suppressing power for wavenumbers k > kj^a- 
Suppression relative to CDM is imprinted on the 
transfer function at matter-radiation equality, so the 
relevant scale is « 9(ma/10“^^ eV) Mpc“^ [3]. 
The relative velocity of DM and baryons also sup¬ 
presses power (in a given patch of size the Vhc coher¬ 
ence length) over a range of wavenumbers centered 
near kybc = aH/vtc |1.3) . The relevant ratio is thus 
kj^al^vbc- At 2 = 40 kybc ~ 300 Mpc“^ and so 


{k‘jl/kybc)U=4o « 0.03(mJ10-22 eV) 


\0.5 


(27) 


I have verified these estimates-by-scale of the rel¬ 
evance of kja versus kybc by numerical solution of 
Eqs. ( 22p^ , with the appropriate average over the 
Vba distribution, following Ref. [13] . The results are 
shown in Fig. at 2 = 40 for various axion masses}^ 


^ I use cosmological parameters of Ref. m .This accounts 
for the difference in normalization to Fig. 2 of Ref. m, 
who used a scale-invariant spectrum and different matter 
content. 


For m < 10“^® eV the axion sound speed cuts off 
power while the effects of vta are only of order a few 
percent. For nia = 10“^® eV the axion sound speed 
cuts off power on scales where Vba is relevant at tens 
of percent in the power. 

Planned 21cm power spectrum experiments may 
be able to detect the Vbc effect for CDM [39MT] . 
The power spectrum shown in Fig. is averaged 
over the vtc distribution, but the power suppression 
caused by relative velocity is highly inhomogeneous. 
The relative velocity modulates the large-scale 21cm 
power spectrum by changing the small-scale matter 
power in regions of order the relative velocity coher¬ 
ence length: the moving background couples small 
and large scales. For example Ref. |39| show that 
this modulation actually increases the 21cm power 
over a range of scales k ^ 0.06 Mpc“^. The increase 
in power compared to the case where relative veloci¬ 
ties are absent is larger than the projected sensitivity 
for telescopes such as the Low Frequency Array. 

As shown in Fig. axion DM with rua ^ 
10“^® eV wipes out small-scale power before mod¬ 
ulation effects due to relative velocity become im¬ 
portant. The small-scale power will be the same in 
different coherence length patches and there will be 
no additional modulation power communicated to 
large scales. For < 10“^® eV the 21cm power on 
large scales will resemble the CDM power in the ab¬ 
sence of relative velocity effects and the power boost 
shown in Ref. (SH) will be absent. A detection of 
the modulation would exclude nia 10“^® eV, an 
improvement of some four orders of magnitude with 
respect to current limits liEI- This prognosis for 
21 cm agrees by order of magnitude with extrapola¬ 
tion from linear theory [42] using the WDM fore¬ 
casts of Ref. [IH] (although these forecasts ignored 
the Vbc effect for WDM, and constraints were driven 
by the reduced collapse fraction, which is expected 
to be similar for axions and WDM jUj). On the 
other hand, the absence of a large-scale modulation 
will imply the existence of some small-scale cut-off 
in the power spectrum at fc < kybc- 

While the small-scale power suppression caused by 
relative velocity is highly inhomogeneous, the power 
suppression caused by the axion sound-speed affects 
all of space. This could lead to further interesting ef¬ 
fects. Finally, and very optimistically, the amplitude 
and shape of the large-scale 21cm power could probe 
axions in the range 10“^® eV < TOq < 10“^® eV, 
where the small-scale power is affected on the same 
scales by both the axion Jeans scale and the Vbc ef¬ 
fect. 
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B. Notes on an SPH treatment 


The identification of the scalar field gradient en¬ 
ergy in the non-relativistic limit with a “quan¬ 
tum force” suggests a new approach to model 
axion/scalar-field DM on small scales. In quantum 
mechanics one can use the hydrodynamic Madelung 
representation to solve the Schrodinger equation by 
introducing the quantum force. This is the so-called 
“synthetic” view of Bohmian mechanics [23]: the 
motion of the hydrodynamic “particles” is not nec¬ 
essarily understood as fundamental, and due to hid¬ 
den variables, although it can nonetheless provide 
insights into processes such as quantum-mechanical 
tunneling 

In axion hydrodynamics, of course, the resem¬ 
blance to Bohmian mechanics is just an analogy. Ax- 
ions are described by a classical field due to the huge 
occupation numbers. One need not worry about 
measurement problems and decoherence. However, 
the hydrodynamic description is one of a condensate: 
although it is a classical held, the dynamics are “very 
quantum” from the particle viewpoint [43] . 

If cosmological structure formation is to be stud¬ 
ied including the effects of the expansion and Hubble 
friction, Eqs. ( Tip^ should be used. If equations 
proposed by Ref. |21j are used, the expansion must 
be correctly accounted for by choice of frame. In the 
absence of the expansion-related terms, one can only 
study already-bound structures. 

It is conceptually simple to consider including the 
axion pressure in SPH. The quantum force in the 
hydrodynamic treatment is simply computed from 
the gradient of Q, and added to the force due to self 
gravity, giving the total force: 


studied, as well as methods for dealing with nodes 
[23] . providing a possible “off-the-shelf cookbook” 
for cosmological simulators. . 

As pointed out by Ref. [2T|, SPH is naturally 
adaptive, smearing the particles over large distances 
in low density regions and providing a natural force 
softening. In the scheme of Ref. [21] the density held 
is represented with smoothing kernels at SPH parti¬ 
cle locations: 

Pa{x) ='^Pa{x^)W{x - Xi',^) , (29) 

i 

Wix; 0 = exp(-|f I V2f), (30) 

(31) 

where ^ is a smoothing parameter. The advantage 
of this method is that density gradients only affect 
the kernel, whose derivative is analytically known. 

In Ref. [ll]’s scheme, ^ is chosen adaptively such 
that there is a hxed mass within the smoothing ker¬ 
nel. This ensures that regions of both high and low 
density are resolved equally well. However, since 
the quantum force can blow up inside low density 
regions such as voids, one might worry about the 
accurate resolution of the velocity field in these re¬ 
gions. Indeed, Ref. HH showed that rapid oscillation 
of the wavefunciton is common everywhere, even in 
voids, implying large fluid velocity. The adaptive 
mesh-refinement in the solution of the Schrodinger 
equation had to be carefully chosen to resolve phase 
oscillations. Similar considerations will apply to op¬ 
timal choices of smoothing in SPH. 

In the continuum limit, the SPH smoothing kernel 
is the position-space Husimi distribution (see, e.g.. 

Ref. [11]) 


F = -a-^V{V + Q). 


(28) 


The gravitational potential, V, is found from the 
density by solving the Poisson equation as usual, 
while Q is computed from gradients of the density 
via Eq. (16). Of course, the correct initial linear 


power spectrum for axion DM must be used. So far, 
so theoretically simple. 

The numerical problem with such a treatment 
arises in two ways. Firstly, Q depends on the Lapla- 
cian of the local density, and derivatives are difficult 
to estimate accurately in multi-scale problems like 
cosmological structure formation. The second diffi¬ 
culty arises because Q can blow up in regions where 
Pa —>■ 0. This will be a problem in particular at in¬ 
terference nodes and in deep voids, where ^ 0. 
In the quantum trajectories literature, many stable 
algorithms for computing the gradients have been 


Ph{x) = J (fx'\ip{x)\'^W{x-x'-,^). (32) 

The smoothing scale, is a free parameter reflecting 
the uncertainty relation. The Husimi distribution 
gives spatial resolution Sx = ^ and momentum res¬ 
olution 6p = (2^)“^. The adaptive SPH smoothing 
described above will make ^ large in regions of low 
density, leading to good velocity resolution in these 
regions and possibly mitigating any issues within 
voids. 

One sees that the choice of SPH smoothing is 
related to choosing the ‘best’ measurements of a 
quantum probability distribution to extract phase- 
space information. Different possible choices of ^ are 
suited to different problems. For example, ^ could 
be chosen using the de Broglie scale, thus increas¬ 
ing resolution in both high density and high velocity 
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regions. Such a smoothing might be useful for study¬ 
ing the phase-space density of axion DM halos |45] 
and the structure of voids in this model. In quantum 
mechanics, the distribution is “nearly classical” if ^ 
can be chosen such that Sx and 6p are each smaller 
than typical structures in the observables. Such a 
criterion could be used to determine when the quan¬ 
tum force is relevant, possibly saving computation 
time in an N-body simulation of this model. 

For further discussion of smoothing issues in the 
Schrodinger picture, see Refs. mUH]; for smoothing 
issues in standard CDM, see e.g. Ref. [iS]. While 
smoothing is ad hoc in CDM, it emerges naturally 
in the scalar held case, along with the adhesion ap¬ 
proximation |25j . 

Since SPH modules are standard in many modern 
N-body codes (as they are used to model the baryon 
sound speed) I hope that this method for modelling 
axion DM on small scales can be implemented rela¬ 
tively easily into cosmological simulations. 


IV. CONCLUSIONS 


I have presented a simple derivation of the non¬ 


linear hydrodynamic equations (NLHEs), Eqs. (14 


16) for axion/scalar held DM derived from a fun¬ 
damental action. They are valid in the large 
occupation-number limit, where axion DM is de¬ 
scribed by the classical held equations. I have taken 
the non-relativistic, Newtonian limit, but have left 
the treatment of the axion energy density and held 
huctuations non-perturbative. The equations pre¬ 
sented differ from the usual Shcrodinger picture for 
DM by the inclusion of expansion and Hubble drag 
effects explicitly, rather than through implicit frame 
choice. The NLHEs are the correct setting to study 
non-relativistic perturbation theory and non-linear 
clustering of axion DM. Two possible applications 


of these equations were outlined. 

The hrst discussed the relative velocity of baryons 
and DM. The effects shown in Fig. will delay for¬ 
mation of hrst stars and thus are expected to affect 
heating of the intergalactic medium and metal en¬ 
richment, and suggests that the 21cm power spec¬ 
trum could tighten constraints on the axion mass 
by several orders of magnitude with respect to the 
current best constraints. Further study of this pos¬ 
sibility is currently underway. If axions play a role 
in the formation of cores in dwarf galaxies, these 
searches will hnd evidence for axions [13 ■ Con¬ 
straining rua ~ 10“^® eV is a theoretically well mo¬ 
tivated goal as this is the “anthropic” boundary for 
axions in string theory |48j . and is also a range 
of masses independently constrained by black hole 
spins [4^[50] . 

The second application discussed an SPH model 
for N-body simulations of axion DM, which is 
an alternative to the adaptive-mesh-rehnement 
approach of Ref. and could be incorporated 

within existing N-body codes. Work on this topic, 
too, is underway. Such simulations will open new 
doors on the study of structure formation with 
axion DM. 
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